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We present a variant of the recently developed quantum corrected model (QCM) for plasmonic 
nanoparticles [Nature Commun. 3, 825 (2012)] using non-local boundary conditions. The QCM 
accounts for electron tunneling in narrow gap regions of coupled metallic nanoparticles, leading 
to the appearance of new charge transfer plasmons. Our approach has the advantages that it 
emphasizes the non-local nature of tunneling and introduces only contact resistance, but not ohmic 
losses through tunneling. Additionally, it can be implemented much easier in boundary element 
method (BEM) approaches. We develop the methodology for the QCM using non-local boundary 
conditions, and present simulation results of our BEM implementation which are in good agreement 
with those of the original QCM. 
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I. INTRODUCTION 

Plasmonics allows to manipulate light at the nanoscale 
and to obtain strong and very confined electromagnetic 
fieidPH. This is achieved by binding light to coher¬ 
ent electron charge oscillations at metal-dielectric inter¬ 
faces, so-called surface plasmons (SPs), sometimes also 
referred to as surface plasmon polaritons. Recent work 
has addressed the question under which conditions a clas¬ 
sical SP description in terms of a local dielectric func¬ 
tion breaks down and quantum-mechanical corrections 
become mandatory. On the one hand, at sharp edges 
and corners of metallic nanoparticles there is a spill-out 
of the electron charge distribution, due to the electron gas 
pressure, which leads to a nonlocal dielectric response 6 9 
causing a blue shift of the SP resonances and a reduction 
of the achievable field enhancements in comparison to lo¬ 
cal description^^. On the other hand, for sub-nanometer 
gaps and sufficiently high field strengths electrons can 
tunnel between neighbor nanoparticles 11 13 leading to 
the emergence of new charge-transfer plasmonJ^. Elec¬ 
tron transfer through larger gaps can occur in molecular 
tunnel junctions 15 . 

From the theoretical side, such quantum corrections 
have been modelled by introducing either modified 
boundary conditions or artificial materials that mimic 
the quantum behaviour. In Ref. 7 the authors showed 
that a non-local dielectric response can be modelled by 
replacing the non-local metal with a composite material, 
comprising a thin dielectric layer on top of a metal with 
local dielectric prope rties. Similarly, in the quantum- 
corrected model 1113 (QCM) an artificial dielectric ma¬ 
terial is filled into the gap region, with a conductivity 
that reproduces the correct tunnel current between two 
neighbour nanoparticles. As the tunnel current typically 
has an exponential dependence with respect to the gap 
distance^, non-planar tunneling gaps must be modelled 
by onion-like shells of materials with different conduc¬ 
tivities. Different materials can be easily introduced in 
volume based simulation approaches, such as finite dif¬ 


ference time domain (FDTD) simulatior l—^ 1 —l 

In this paper we show how to simulate tunneling effects 
within a boundary element method (BEM) approach 19 21 
by introducing modified non-local boundary conditions. 
While the consideration of additional materials is com¬ 
putationally cheap in volume based simulations, it be¬ 
comes computationally very demanding in BEM simula¬ 
tions, since usually a large number of different material 
layers is needed to resolve the exponential tunnel cur¬ 
rent dependence. In contrast, the consideration of mod¬ 
ified boundary conditions in a QCM variant has virtu¬ 
ally no impact on the performance of BEM simulations 
compared to conventional ones. We will show that both 
approaches, either the consideration of artificial materi¬ 
als or modified non-local boundary conditions, give sim¬ 
ilar results. From a conceptual point of view, non-local 
boundary conditions have the advantage that they em¬ 
phasize the non-local behaviour of the tunneling process 
and tunnel currents do not suffer from ohmic losses but 
are only governed by contact resistance, a finding known 
for a long time in the field of mesoscopic electron trans¬ 
port!^. 

II. THEORY 

Fi gure 1 (a) shows the basic principle of the original 
QCMUDSI (i n the following denoted as volume QCM) at 
the example of two nanoparticles separated by a small 
gap of sub-nanometer size. When an electric field E is 
applied across the gap, a tunnel current 

Jt = &t E (1) 

starts to flow, where a t is the tunnel conductivity that 
can be either obtained from first principles or effec¬ 
tive modelcalculations of various degrees of sophistica¬ 
tion 11 13 23 24 . To mimic such tunnel currents, within the 
quantum corrected model one introduces in the gap re¬ 
gion an effective, homogeneous medium e^t with a con¬ 
ductivity chosen to yield the correct tunnel current (we 









2 




FIG. 1: (Color online) Schematics of the quantum corrected 
mod el (QC M). (a) Volume based implementation of Esteban 
et al! 11 13 where artificial dielectric materials are placed in¬ 
side the gap. The conductivities of these materials are set to 
the gap-size dependent tunnel conductivities, (b) Boundary 
element based implementation of this work, with non-local 
artificial boundary conditions which are chosen in order to 
obtain the proper tunnel current between boundary positions 
s a and Sb- The inset indicates the pillbox (with outer surface 
normal n a ) over which Gauss’ law is integrated to obtain the 
artificial boundary conditions. For details see text. 

adopt the notation of Ref. [19] and denote the dielectric 
functions in- and outside the nanoparticle with £i and £ 2 , 
respectively). This approach has a number of advantages: 
first, it can be easily implemented in volume based sim¬ 
ulation approaches, such as FDTD; second, the descrip¬ 
tion in terms of a local current distribution guarantees 
that charge is conserved, i.e., the charge that leaves one 
nanoparticle must be transferred via the junction to the 
other nanoparticle. On the other hand, the approach has 
a number of conceptual difficulties: the current is sub¬ 
ject to ohmic losses, contrary to the purely contact-like 
resistivity of quantum tunneling; additionally, current is 
not only induced by electric fields parallel the nanoparti¬ 
cle connection, such as one would expect for tunnel cur¬ 
rents, but also by perpendicular fields. In most cases of 
interest these are no serious shortcomings, since fields in 
gap regions practically always point along the nanopar¬ 
ticle connection, and the tunnel junction is typically so 
narrow that ohmic losses are of only minor importance. 

We will next rephrase the QCM in terms of modified 
boundary conditions which are much better suited for 
BEM implementations. Our starting point is Gauss’ law 
integrated over the small pillbox indicated in Fig. 1(b), 

V • D dr = j) D • da = 47 r J pdr 

4tt r zi-7 to r 

= — / V J t dr =- <bJ t -da, ( 2 ) 

iuj J uj J 


Here a and b denote the left and right nanoparticle, re¬ 
spectively. The last term accounts for the charge trans¬ 
ferred from position s a to s 5 through quantum tunneling 
(i.e., the loss or gain of charge in the pillbox over which 
Gauss’ law is integrated). Similarly to Eq. ([T]) we as¬ 
sume that the current is proportional to the tunnel con¬ 
ductivity a t and the average of the electric field along 
the outer surface normal directions n a ^ [as n a and n& 
in the gap region are approximately antiparallel, E in 
Eq. © receives a negative sign]. Note that this choice 
is by no means unique. We could alternatively assume 
J at = cF t {E 2a + E 2 b )/2 or Jat = cr t E[(s a + s b )/ 2 ], In 
all cases charge remains conserved since the current J at 
leaving particle a at position s a is always the opposite 
to the current J^t entering particle b , and vice versa. 
However, the consideration of solely normal currents 
has the advantage that only the boundary condition of 
the dielectric displacement needs to be modified, whereas 
the boundary condition for the parallel magnetic field re¬ 
mains unaltered because of our neglect of parallel tunnel 
currents. 

Eq. © is the central result of this work. It replaces the 
consideration of artificial dielectric materials through an 
artificial bou ndary condition. Contrary to the QCM of 
Esteban et al ! n * 13 ( our approach describes quantum tun¬ 
nel as a genuine non-local process and thus does not suffer 
from ohmic losses in the tunnel junction. It can be also 
easily extended to molecular tunnel junction by lumping 
all microscopic details about the microscopic tunneling 
process into an effective a t value. As regarding the role 
of normal and parallel electric fields in tunneling, both 
models are comparably arbitrary but could be further 
refined. However, since in narrow gap regions the plas- 
monic nearfields preferentially point along the interpar¬ 
ticle connection, the detailed E 1 - and E^ behavior of cr t 
is usually completely irrelevant. 

In Appendix [A] we show how to modify the BEM ap¬ 
proach of Ref. [19] to account for quantum tunneling, and 
present the working equations that can be implemented 
within the MNPBEM toolboxPlX] 


III. RESULTS 

We start by considering in accordance to Refs. II 11131 
the case of two spheres with a gap in the sub-nanometer 
regime. For the dielectric function we take a Drude-type 
form e((jj) = £ 0 — uj 2 /(uj 2 + iuj'j) for gold, £2 = 1 for the 
embedding medium, and 


where dr and da denote volume and surface integrations, 
respectively, and we have used the Fourier transformed 
continuity equation to relate p t to J t (we use Gaussian 
units throughout). We now make the following ad-hoc as¬ 
sumption for the boundary condition of the normal com¬ 
ponent of the dielectric displacement 
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£2 t{£) — 1 + 


47 r ia t (£) 


cr t (£) = Jm 


UJZ. 


uj 2 + iujj p e^ 


for the tunnel material 13 . Here £0 = 10, uj p = 9.065 eV, 
7 P = 0.0708 eV, and £ c = 0.04 nm, and we consider only 
purely imaginary conductivity corrections for the tunnel 
material. These model parameters provide a good fit to 
experimental datsP 5 for photon energies below 2 eV but 
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FIG. 2: (Color online) Comparison of volume quantum cor¬ 
rected model (QCM) of Esteban et al! n | 13 1 with the bound¬ 
ary QCM of this work. We use two spheres with diameters 
of 50 nm and a Drude-type dielectric function representative 
of gold, and a single layer of artificial tunnel material. The 
light polarization is along the nanoparticle connection. The 
material covers a distance range between the gap size d gap 
and dgap + 0.2 nm, and the artificial dielectric function is 
£ 2 t(d gap + 0.1 nm). The figure shows the gap-size dependent 
extinction cross section (offset for clarity, gap distance given 
on left axis) for the volume QCM and compares them with 
results of the boundary QCM. In the latter approach, we con¬ 
sider quantum tunneling in the same distance window as in 
the volume QCM, and set the tunneling dielectric function to 
the same value as in the volume QCM. 


underestimate dielectric losses above 2 eV where d-band 
scatterings set in. Nevertheless, in this work we keep 
the Drude description to facilitate the comparison with 
Refs. 1111131 The frequency dependence and details of a t 
are subject of ongoing research efforts, the parametriza- 
tion of Eq. @ has been motivated by static tunneling 
calculations including image charge effects as well as by 
time-dependent density functional theory calculations^!, 
related work has employed theory developed for optical- 
assisted tunneling in the microwave domain 23 or dia¬ 
grammatic expansions for the ac conductance through 
inclusion of higher-order electron-plasmon interactions 24 . 
As the primary goal of this work is the derivation and 
implementation of a boundary QCM using a suitable a t 
parametrization, we will here not further elaborate on 
this point. 

Fig. 2 compares for a single artificial tunnel material in 
between the two spheres (see inset) the extinction cross 
sections for different gap distances d gap . The material 
covers the distance range from d gap to d gap + 0.2 nm and 
the dielectric function £ 2 t{dg&p + 0.1 nm) is evaluated at 
the average distance. For the boundary QCM we use 
the same value for £21 and connect boundary elements 
of the two neighbour spheres within the same distance 
range 26 . With this, we are able to compare the volume 


FIG. 3: (Color online) Volume and boundary QCM for the 
same spheres as in Fig. 2 and for d gap = 0.075 nm. In the vol¬ 
ume QCM we consider an onion-like sequence of five materials 
e(£), with i covering the region from d gap to d gap + 0.4 nm. 
In the boundary QCM we use the e{tj values for the respec¬ 
tive boundary element distances. Volume QCM1 refers to the 
model of Ref. [lT| and volume QCM2 to a simulation where 
the light excitation and the scattered far fields are computed 
without the artificial materials. Boundary QCM1 refers to 
simulations where opposite boundary elements of the flipped 
spheres are connected (with a refined mesh at the poles), and 
boundary QCM2 to a simulation where the respective closest 
boundary elements of the neighbour spheres are connected. 


and boundary QCM directly. As can be seen in the fig¬ 
ure, both volume and boundary QCM give practically 
identical results over the entire range of gap distances 
where tunneling sets in. Tunneling is evidenced by the 
disappearance of the lowest plasmon peak around 1.8 eV 
with decreasing gap distance, and the onset of the charge 
transfer peak around 0.8 eV. Similarly to the extinction 
spectra, also the field enhancements in the gap region 
(not shown) computed within the volume and boundary 
QCM are in almost perfect agreement. It is gratifying to 
see that the volume and boundary QCM models compare 
so well. 

Next, we show in Fig. 3 results for the full QCM simula¬ 
tions for the same setup as in Fig. 2 and for d gap = 0.075 
nm. For the volume QCM we use five layers of artifi¬ 
cial materials, covering the distance range from d gap to 
d gap + 0.2 nm, and for the boundary QCM we use for 
£ 2 t(f) the respective distances i between opposite bound¬ 
ary elements. Note that we use for both spheres the same 
boundary meshes with a refined discretization at one of 
the poles®! and simply flip and displace the spheres to 
obtain the dimer structure shown in the inset. Again we 
find good agreement between the volume and boundary 
QCM, although the volume QCM leads to a more pro¬ 
nounced exctinction peak of the charge transfer plasmon. 

We believe that this is an artefact caused by our BEM 
implementation of the volume QCM. The BEM approach 
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FIG. 4: (Color online) Extinction cross section for a dimer, us¬ 
ing a classical electrodynamic (gray, dashed line) and a QCM 
simulation (blue line) with polarization along the nanoparti¬ 
cle connection, as well as a QCM simulation for a trimer (red 
line). The sphere diameters are 50 nm and the gap distances 
are 0.1 nm. For the trimer, the optical spectra do not de¬ 
pend on the polarization direction of the incoming light (light 
propagation direction perpendicular to trimer plane). 


of Garcia de Abajo and Howie matches electromag¬ 
netic potentials at mater ial boundaries in order to solve 
Maxwell’s equations 19 20 . In this approach, an external 
plane wave excitation only excites materials connected 
with the embedding medium (in the gap region the out¬ 
ermost material is the last layer of artificial tunneling ma¬ 
terial) and the excitation is then passed to the inner lay¬ 
ers through the solution of Maxwell’s equations 19 . While 
this causes typically no problems, it becomes computa¬ 
tionally demanding for the inhomogeneous tunnel mate¬ 
rial which is modelled through closely spaced onion-like 
layers. In our simulations we had problems to get fully 
converged results when increasing the number of layers, 
probably due to artificial reflections and transmissions of 
the incoming light at the layer interfaces. When we con¬ 
sider the tunneling materials only in the BEM solutions 
and (artificially) neglect them in the light excitation (see 
simulation results with diamond symbols) we obtain for 
the charge transfer peak perfect agreement between vol¬ 
ume and boundary QCM. Also the (minor) differences 
at higher energies are probably due to implementation 
problems of the volume QMC within the BEM approach. 

The squares in Fig. 3 report results of a slight variant 
of the boundary QCM. Here we do not connect oppo¬ 
site boundary elements (as one can only do for flipped 
nanoparticles) but connect the closest boundary elements 
of the two nanoparticles. Apparently, such an approach 
also works for nanoparticle arrangements with a lower 
degree of symmetry. As one infers from a comparison of 
the boundary QCM1 and QCM2 results, these two ap¬ 
proaches are in perfect agreement. 


As a final example, in Fig. 4 we show results for a 
symmetric trimer structure consisting of three spheres, 
demonstrating that simulations of more complicated 
nanoparticles and nanoparticle arrangements can be eas¬ 
ily performed with our BEM approach. For the trimer 
structure we again observe the appearance of the charge 
transfer plasmon peak. Due to the triangular symme¬ 
try, the extinction cross sections do not depend on the 
polarization of the incoming light (propagating perpen¬ 
dicularly to the trimer plane). 

IV. SUMMARY AND CONCLUSIONS 

To summarize, we have presented a variant of the 
quantum corrected model (QCM) where tunneling is ac¬ 
counted for by the consideration of non-local boundary 
conditions. This approach has the advantage that it em¬ 
phasizes the non-local nature of tunneling and does not 
introduce artificial ohmic tunnel losses. We have devel¬ 
oped the methodology for implementing the boundary 
QCM within a boundary element method (BEM) ap¬ 
proach, and have presented simulation results which have 
compared well with results of the original volume QCM. 
Minor differences between the two approaches have been 
attributed to intrinsic difficulties of our BEM scheme to 
properly implement a volume QCM. We believe that the 
volume and boundary QCM are closely related, but the 
availability of a different approach might be beneficial for 
conceptual reasons as well as for BEM implementations. 

Our approach might prove particularly useful for 
molecular tunnel junctions with larger gap sizes. Also 
supplementing the QCM through inclusion of non-local 
effects in the dielectric metal function, through modi¬ 
fied boundary conditions, should be relatively straight¬ 
forward. Future work will also address the possibilities to 
compute the tunnel conductivities through ab-initio cal¬ 
culations and to submit the pertinent tunnel parameters 
to classical electrodynamic simulations including quan¬ 
tum corrections. 
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Appendix A 

Here we show how to implement the non-local quan¬ 
tum tunneling of Eq. |3| in the BEM approach of Garcia 
de Abajo and Howie®(in the following we refer to the 
equations of this work with a preceding G). Importantly, 
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we can carry over most results with the only exception 
of Eqs. (G17,G18) which become modified through the 
nonlocal boundary condition. 

The continuity of the scalar and vector potentials 0 
and A read [Eqs. (G 10 ,G 11 )] 

G\(J\ - G 2 (J 2 = 02 - 01 = F 
G 1 h 1 -G 2 h 2 = A\ — A\ = a , 

where G i and G 2 denote the Green functions inside and 
outside the nanoparticle, and a and h are artificial sur¬ 
face and current distributions at the particle boundary 
which are chosen such that the boundary conditions of 
Maxwell’s equations are fulfilled. 0 e and A e are the 
scalar and vector potentials of an external excitation, 
such as a plane wave. For further details see Refs. [191201 
The continuity of the magnetic field becomes [see also 
Eq. (G14)] 

Hihi — H 2 h 2 — ik ft ( e\G\(J\ — e 2 G 2 a 2 ) = ol 

with H\ 2 being the surface derivative of Gy 2 taken at 
the particle in- or outside, and a' is defined through 
Eq. (G15). For the continuity of the normal dielectric 
displacement we get 

EiHicji- e 2t H 2 a 2 -ik {sih • Gihi - e 2t h • G 2 h 2 ) = D ef , 


Yet, the derivation of the BEM equations is not too dif¬ 
ferent. 

First, we use 

G\(j\ = G 2 (j 2 + ip 
G\h\ = G 2 h 2 + a 

to replace in the continuity equation (G14) of the mag¬ 
netic field <ti, hi by a 2 , h 2 , 

(Ei - E 2 ) G 2 h 2 -ikh (s 1 - e 2 ) G 2 a 2 = a , 

with Ei = HiG ^ 1 , E 2 = H 2 G 2 X and a = a! — Eia + 
ik heiip. The continuity of the normal dielectric displace¬ 
ment becomes 

(^lEi — s 2 t^ 2 ) G 2 (j 2 — ik (^i — e 2 t) ft • G 2 h 2 = D e , 

with D e = D e ' — £iEi<p + iksih • a. We can use the 
continuity equation for the magnetic field to express the 
surface current h 2 in terms of a 2 , 

G 2 h 2 = A -1 [ik h(e i - e 2 )G 2 a 2 + a] , (Al) 
with A = Ei — E 2 . Inserting this expression into the con¬ 
tinuity equation for the normal dielectric displacement we 
finally obtain 


with 


D 


ef 


£i {ik h • A\ — 0®') — s 2t {ik h ■ A\ — 0^) . 


6iEi - 5 2t E 2 + k 2 {e i - e 2t )h • A 1 h{s 1 - e 2 ) 
= D e - 1 - ik{e i — e 2t )h • A _1 a . 


G 2 cf 2 

(A2) 


Here 0f 2 ' denote the surface derivatives of the external 
scalar potentials, and s 2t = e 2 + (4 Tti(j t /uj) is a non-local 
dielectric function accounting for quantum tunneling, see 
Eq. <©• Because s 2t is nonlocal and connects points s a 
and St, through tunneling, it cannot be commuted with 
the Green functions as in the original BEM approach 19 . 


Equations ( |A 1 | ) and ( |A 2 | ) are the two working equa¬ 
tions of our BEM approach which can be solved through 
matrix inversion. Once the surface charges and currents 
cr 2 and h 2 are known for a given external excitation, one 
can compute the electrodynamic potentials and fields ev¬ 
erywhere else. 
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